function [theta_est,covar,p_unob_het] = AM_Estimation2(FID_Grid,States,AgeState,Action,pproc,opt)
% Estimate dynamic model using the NPL method proposed by Aguirregabiria
% and Mira (2003). 
% Error correction for underobserved replanting based on Arcidiacono and
% Miller (2011).

% Before beginning the estimation per se, we need the matrix with exogenous 
% transition probabilities
Regions = unique(States.RegionsForExoStates);
Nr = length(Regions);
NumR = length(Regions);
for k = 1:Nr % State transition:
    reg = Regions(k);
    [F{reg}] = StateTrans(States.AllStatesMap,pproc.P{reg});
%      [F] = StateTrans(States.AllStatesMap,pproc.P{reg});
%      save([opt.ScratchFolder 'F' num2str(reg) '.mat'],'F','-v7.3');
%      clear F
end

% Prepare data for estimation. In the likelihood function, we will have observations stacked 
% (and not each year in a single column)
TAll = [2003:2013]; % All years used.
TUse = opt.Years;
TNotUse = setdiff(TAll,TUse);
StackedAction = []; StackedSingleState = []; StackedFieldID = [];
for t=1:length(TAll)
    year = TAll(t);
    if ~any(TNotUse == year) % If year is not in the NotUse set, that is, if the year is to be used... 
        StackedAction        = [StackedAction; Action(:,t)];
        StackedSingleState   = [StackedSingleState; States.SingleStateIdx(:,t)];
        StackedFieldID       = [StackedFieldID; FID_Grid];
    end
end

% field ID in indexes for estimation:
[~, ~, IndFieldID] = unique(StackedFieldID);

% Set up initial parameter vector:
Nf     = size(States.ExoFixedStates,2)-3; % # exo cost variables
theta0 = [opt.ini.gamma;opt.ini.thetaR;opt.ini.PsiE*ones(5,1);opt.ini.PsiR;opt.ini.PsiA;opt.ini.kappa; ...
    ones(2*Nf,1)*opt.ini.beta];

options = optimoptions('fminunc',                                        ...
                       'Algorithm',                 'quasi-newton',      ...        
                       'GradObj',                   'on',                ...
                       'DerivativeCheck',           'off',               ...
                       'Display',                   'off',               ...
                       'Diagnostics',               'off',               ...
                       'TolFun',                     1e-12,              ...
			           'TolX',                       1e-10,              ...
                       'MaxIter',                    700);

tol = opt.NPLtol;                   
% Start NPL k-steps iteration:
disp(' ')
if opt.IsBoot == 0
    ES = exist([opt.SaveResultsAs '_step.mat']);
    if opt.Restart == 1 && ES == 2
        load([opt.SaveResultsAs '_step.mat']); % This should have pOld, thetaOld and iter.
        disp(' Previous result found and loaded ... ')
    else
        % Run first step to get initial CCPs:
        disp(' No previous result found, starting from scratch ... ')
        disp(' ')
        disp(' Running first stage ... ')
        [pOld,theta_fs] = FirstStage(Action,AgeState,States,pproc);
        p_uh_old = 0.85;
        thetaOld = theta0;
        iter = 0;
    end
else 
    iter = 0;
    % Bootstrap estimation
    % Should have thetaOld, pOld loaded in opt structure.
    thetaOld = opt.thetaOld;
    pOld     = opt.pOld;
    p_uh_old = opt.p_uh_old;
    tol = 10^-4; % Reduce tol for bootstrap for speed (decrease precision: > variance)
end


disp(' ')
err = 1; %max_iter = 1000;

fprintf(1,'%20s',  'Exit (fminunc)'); fprintf(1,'%15s',  'Err'); fprintf(1,'%20s',  'Time Elapsed'); fprintf(1,'%10s',  'Iter'); fprintf(1,'%13s', 'p_uh'); fprintf(1,'\n');

while err > tol && iter < opt.max_iter
    
    
    tic
    
    [UHprob,p_uh_new] = UpdateUHProb2(p_uh_old,States,AgeState,Action,pOld);
    
    StackedUHprob = repmat(UHprob,[length(TUse) 1]);
    
    [thetaNew,~,exit] = fminunc(@(x) LikAM(x,F,States,StackedUHprob,StackedAction,StackedSingleState,IndFieldID,pproc,pOld,opt,0),thetaOld,options);

    % Get new pNew:
    for r = Regions'
        reg = r;
        [U{r},dUdTheta{r}] = Payoff(States.AllStatesMap,States.MidPointsFixedStates,States.FixedStatesMap,States.SRDMidPoints,thetaNew,pproc.y{reg});
        % Renewal choice we will work with is choice 1 (Start sugar or Replant)
        
        %load([opt.ScratchFolder 'F' num2str(r) '.mat'],'F');
        % Compute Log(P): rows -> states, column -> Action.
        for a=1:3
            %ExpA(:,a) = exp(U{r}(:,a) - opt.rho*F{r}{a}*log(pOld{r}(:,2)));
            ExpA(:,a) = exp(U{r}(:,a) + opt.rho*F{r}{a}*( U{r}(:,2) - log(pOld{r}(:,2))) );
        end
        %clear F
        
        DenA = sum(ExpA,2);
        pNew{r} = bsxfun(@ldivide,DenA,ExpA);
        
    end
    
    % Calculate error:
    errTheta = max(abs(thetaNew-thetaOld));
    r_idx = 0;
    for r = Regions'
        r_idx = r_idx + 1;
        errP(r_idx) = max(max(abs(pOld{r}-pNew{r}))); 
    end
    err = max([errTheta,errP]);
    
    TimeElap = toc;

    fprintf(1,'%20.0f', exit);
    fprintf(1,'%15.6f', err);
    fprintf(1,'%20.6f', TimeElap);
    fprintf(1,'%10.0f', iter);
    fprintf(1,'%10.3f', p_uh_old);
    fprintf(1,'\n');
    
    
    % Display convergence:
    disp('  ')
    
    iter = iter + 1;
    pOld = pNew;
    thetaOld = thetaNew;
    p_uh_old = p_uh_new;
    
    if opt.Restart == 1 && opt.IsBoot == 0
        save([opt.SaveResultsAs '_step.mat'],'thetaOld','pOld','p_uh_old','iter');
    end

end

theta_est  = thetaNew;
p_unob_het = p_uh_new;

% Compute covariance matrix:
[~,~,covar] =  LikNPL(theta_est,F,States,StackedAction,StackedSingleState,IndFieldID,pproc,pOld,opt,1);


end